namespace Eigen {

namespace internal {

template<typename Scalar>
void
r1updt(Matrix<Scalar, Dynamic, Dynamic>& s,
	   const Matrix<Scalar, Dynamic, 1>& u,
	   std::vector<JacobiRotation<Scalar>>& v_givens,
	   std::vector<JacobiRotation<Scalar>>& w_givens,
	   Matrix<Scalar, Dynamic, 1>& v,
	   Matrix<Scalar, Dynamic, 1>& w,
	   bool* sing)
{
	typedef DenseIndex Index;
	const JacobiRotation<Scalar> IdentityRotation = JacobiRotation<Scalar>(1, 0);

	/* Local variables */
	const Index m = s.rows();
	const Index n = s.cols();
	Index i, j = 1;
	Scalar temp;
	JacobiRotation<Scalar> givens;

	// r1updt had a broader usecase, but we don't use it here. And, more
	// importantly, we can not test it.
	eigen_assert(m == n);
	eigen_assert(u.size() == m);
	eigen_assert(v.size() == n);
	eigen_assert(w.size() == n);

	/* move the nontrivial part of the last column of s into w. */
	w[n - 1] = s(n - 1, n - 1);

	/* rotate the vector v into a multiple of the n-th unit vector */
	/* in such a way that a spike is introduced into w. */
	for (j = n - 2; j >= 0; --j) {
		w[j] = 0.;
		if (v[j] != 0.) {
			/* determine a givens rotation which eliminates the */
			/* j-th element of v. */
			givens.makeGivens(-v[n - 1], v[j]);

			/* apply the transformation to v and store the information */
			/* necessary to recover the givens rotation. */
			v[n - 1] = givens.s() * v[j] + givens.c() * v[n - 1];
			v_givens[j] = givens;

			/* apply the transformation to s and extend the spike in w. */
			for (i = j; i < m; ++i) {
				temp = givens.c() * s(j, i) - givens.s() * w[i];
				w[i] = givens.s() * s(j, i) + givens.c() * w[i];
				s(j, i) = temp;
			}
		} else
			v_givens[j] = IdentityRotation;
	}

	/* add the spike from the rank 1 update to w. */
	w += v[n - 1] * u;

	/* eliminate the spike. */
	*sing = false;
	for (j = 0; j < n - 1; ++j) {
		if (w[j] != 0.) {
			/* determine a givens rotation which eliminates the */
			/* j-th element of the spike. */
			givens.makeGivens(-s(j, j), w[j]);

			/* apply the transformation to s and reduce the spike in w. */
			for (i = j; i < m; ++i) {
				temp = givens.c() * s(j, i) + givens.s() * w[i];
				w[i] = -givens.s() * s(j, i) + givens.c() * w[i];
				s(j, i) = temp;
			}

			/* store the information necessary to recover the */
			/* givens rotation. */
			w_givens[j] = givens;
		} else
			v_givens[j] = IdentityRotation;

		/* test for zero diagonal elements in the output s. */
		if (s(j, j) == 0.) {
			*sing = true;
		}
	}
	/* move w back into the last column of the output s. */
	s(n - 1, n - 1) = w[n - 1];

	if (s(j, j) == 0.) {
		*sing = true;
	}
	return;
}

} // end namespace internal

} // end namespace Eigen
